Full scRNAseq pipeline

Author

Pierre-Luc Germain & Anthony Sonrel

Loading necessary libraries

suppressPackageStartupMessages({
  library(SingleCellExperiment)
  library(scran) # sc analysis
  library(scater) # sc QC and plotting
  library(batchelor) # batch correction methods
  library(scDblFinder) # doublet detection
  # library(sctransform) # variance-stabilizing transformation (optional)
  library(muscat) # differential expression analysis
  library(BiocParallel) # handling multi-threading
  library(BiocNeighbors) # specifying params for kNN
  library(igraph) # for manual graph clustering
  library(sechm) # plotting
  library(ggplot2) # plotting
  library(patchwork) # to combine
  library(grid)
})
Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
glmmTMB was built with TMB version 1.9.6
Current TMB version is 1.9.10
Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)

Preprocessing & Clustering

Loading the data

## if you haven't downloaded it yet:
download.file("http://imlspenticton.uzh.ch/teaching/STA426/week13.SCE.rds", dest="week13.SCE.rds")
# note that by default downloads in R timeout after 1min, so
# when downloading large objects you might need to increase this limit:
# options(timeout=1200) # timeout after 20min, instead of 1min...
sce <- readRDS("week13.SCE.rds")

# Have a look at the object
sce
class: SingleCellExperiment 
dim: 27998 12147 
metadata(0):
assays(1): counts
rownames(27998): ENSMUSG00000051951.Xkr4 ENSMUSG00000089699.Gm1992 ...
  ENSMUSG00000096730.Vmn2r122 ENSMUSG00000095742.CAAA01147332.1
rowData names(2): ENSEMBL SYMBOL
colnames(12147): AAACCTGCAAAGCGGT-1.LC017 AAACCTGCACGGACAA-1.LC017 ...
  TTTGTCACATCTATGG-1.LC025 TTTGTCAGTCATATGC-1.LC025
colData names(3): sample_id barcode group_id
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):
# As any matrix it has dimensions
dim(sce)
[1] 27998 12147
# To get the count assay:
counts(sce)[41:45,45:50] #same as assays(sce)$counts
5 x 6 sparse Matrix of class "dgCMatrix"
                                 AAGGAGCAGGAATTAC-1.LC017
ENSMUSG00000042501.Cpa6                                 .
ENSMUSG00000048960.Prex2                                2
ENSMUSG00000057715.A830018L16Rik                        .
ENSMUSG00000097171.Gm17644                              .
ENSMUSG00000101314.Gm29663                              .
                                 AAGGAGCGTAGGAGTC-1.LC017
ENSMUSG00000042501.Cpa6                                 .
ENSMUSG00000048960.Prex2                                .
ENSMUSG00000057715.A830018L16Rik                        2
ENSMUSG00000097171.Gm17644                              .
ENSMUSG00000101314.Gm29663                              .
                                 AAGGCAGTCTCTGAGA-1.LC017
ENSMUSG00000042501.Cpa6                                 .
ENSMUSG00000048960.Prex2                                .
ENSMUSG00000057715.A830018L16Rik                        1
ENSMUSG00000097171.Gm17644                              .
ENSMUSG00000101314.Gm29663                              .
                                 AAGGTTCAGCTAGTTC-1.LC017
ENSMUSG00000042501.Cpa6                                 .
ENSMUSG00000048960.Prex2                                .
ENSMUSG00000057715.A830018L16Rik                        2
ENSMUSG00000097171.Gm17644                              .
ENSMUSG00000101314.Gm29663                              .
                                 AAGGTTCAGCTCTCGG-1.LC017
ENSMUSG00000042501.Cpa6                                 1
ENSMUSG00000048960.Prex2                                .
ENSMUSG00000057715.A830018L16Rik                        2
ENSMUSG00000097171.Gm17644                              .
ENSMUSG00000101314.Gm29663                              .
                                 AAGTCTGGTGTTTGTG-1.LC017
ENSMUSG00000042501.Cpa6                                 1
ENSMUSG00000048960.Prex2                                .
ENSMUSG00000057715.A830018L16Rik                        2
ENSMUSG00000097171.Gm17644                              .
ENSMUSG00000101314.Gm29663                              .
# To get the cell names and cell metadata
head(colnames(sce))
[1] "AAACCTGCAAAGCGGT-1.LC017" "AAACCTGCACGGACAA-1.LC017"
[3] "AAACCTGTCGCAAACT-1.LC017" "AAACGGGGTATTCGTG-1.LC017"
[5] "AAACGGGGTCGCGGTT-1.LC017" "AAACGGGGTGGTTTCA-1.LC017"
head(colData(sce))
DataFrame with 6 rows and 3 columns
                         sample_id            barcode group_id
                          <factor>        <character> <factor>
AAACCTGCAAAGCGGT-1.LC017 LC017_LPS AAACCTGCAAAGCGGT-1      LPS
AAACCTGCACGGACAA-1.LC017 LC017_LPS AAACCTGCACGGACAA-1      LPS
AAACCTGTCGCAAACT-1.LC017 LC017_LPS AAACCTGTCGCAAACT-1      LPS
AAACGGGGTATTCGTG-1.LC017 LC017_LPS AAACGGGGTATTCGTG-1      LPS
AAACGGGGTCGCGGTT-1.LC017 LC017_LPS AAACGGGGTCGCGGTT-1      LPS
AAACGGGGTGGTTTCA-1.LC017 LC017_LPS AAACGGGGTGGTTTCA-1      LPS
# To get the gene names and genes metadata
head(rownames(sce))
[1] "ENSMUSG00000051951.Xkr4"    "ENSMUSG00000089699.Gm1992" 
[3] "ENSMUSG00000102343.Gm37381" "ENSMUSG00000025900.Rp1"    
[5] "ENSMUSG00000109048.Rp1"     "ENSMUSG00000025902.Sox17"  
head(rowData(sce))
DataFrame with 6 rows and 2 columns
                                      ENSEMBL      SYMBOL
                                  <character> <character>
ENSMUSG00000051951.Xkr4    ENSMUSG00000051951        Xkr4
ENSMUSG00000089699.Gm1992  ENSMUSG00000089699      Gm1992
ENSMUSG00000102343.Gm37381 ENSMUSG00000102343     Gm37381
ENSMUSG00000025900.Rp1     ENSMUSG00000025900         Rp1
ENSMUSG00000109048.Rp1     ENSMUSG00000109048         Rp1
ENSMUSG00000025902.Sox17   ENSMUSG00000025902       Sox17
table(sce$sample_id, sce$group_id)
           
              WT  LPS
  LC016_WT  1730    0
  LC017_LPS    0 1027
  LC019_WT  1064    0
  LC020_LPS    0 1468
  LC022_WT  1806    0
  LC023_LPS    0 1546
  LC025_WT  1259    0
  LC026_LPS    0 2247

More about SCE class

QC

# scater / scuttle
# get mitochondrial genes:
mito <- grep("mt-", rownames(sce), value = TRUE)
head(mito)
[1] "ENSMUSG00000064341.mt-Nd1"  "ENSMUSG00000064345.mt-Nd2" 
[3] "ENSMUSG00000064351.mt-Co1"  "ENSMUSG00000064354.mt-Co2" 
[5] "ENSMUSG00000064356.mt-Atp8" "ENSMUSG00000064357.mt-Atp6"
# get QC metrics:
sce <- addPerCellQC(sce, subsets=list(Mt=mito), percent.top=c(5,10))
sce <- addPerFeatureQC(sce)

# we plot some of the metrics
qc <- as.data.frame(colData(sce))
ggplot(qc, aes(subsets_Mt_percent)) + geom_histogram() + facet_wrap(~sample_id)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(qc, aes(log10(sum))) + geom_histogram() + facet_wrap(~sample_id)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(qc, aes(log10(sum), log10(detected))) + geom_point() + geom_density2d()

# we set thresholds on the library sizes and detection rate:
sce$qc.out <- isOutlier(log(sce$sum),nmads=3,batch=sce$sample_id) |
              isOutlier(log(sce$detected),nmads=3,batch=sce$sample_id)
# a fancier job would be accomplished by https://github.com/wmacnair/SampleQC
table(sce$qc.out)

FALSE  TRUE 
12115    32 
# get rid of seldom detected genes
sce <- sce[rowData(sce)$detected>=4,]

# we flag doublets: doublets are droplets that capture more than one cell and this might be a problem when it captures more than one cell-type
sce <- scDblFinder(sce, samples="sample_id", BPPARAM=MulticoreParam(6))
table(sce$scDblFinder.class)

singlet doublet 
  11484     663 

Normalization & reduction

# (fast) standard log-normalization
sce <- logNormCounts(sce)

# get high-variable genes (normally we'd take 3000, here we'll go skim for speed)
hvg <- getTopHVGs(sce, n=1000)

# alternative: variance-stabilizing transformation
# vst <- suppressWarnings(sctransform::vst(counts(sce)))
# logcounts(sce) <- vst$y
# # get highly-variable genes
# hvg <- row.names(sce)[order(vst$gene_attr$residual_variance,
#                             decreasing=TRUE)[1:1000]]


# run PCA
sce <- runPCA(sce, ncomponents=50, subset_row=hvg)

#How many componentes do you take? before the used the elbo, but it is not a good way to proceed. 

# check the variance explained by the PCs:
pc.var <- attr(reducedDim(sce),"percentVar")
plot(pc.var, xlab="PCs", ylab="% variance explained")

# restrict to the first 20 components:
reducedDim(sce) <- reducedDim(sce)[,1:20]

# run and plot 2d projections based on the PCA
# to improve performance, use Annoy kNN approximation:
sce <- runTSNE(sce, dimred="PCA", BNPARAM=AnnoyParam())
sce <- runUMAP(sce, dimred="PCA", n_neighbors=30, BNPARAM=AnnoyParam())

# you can compare the 2d embeddings:
# sleepwalk::sleepwalk( as.list(reducedDims(sce)[c("TSNE","UMAP")]),
#                       featureMatrices=reducedDims(sce)[["PCA"]] )

#it is an impossible task to do not introduce distorsions in these methods of dimensionality reduction.  Because to perserve the neighbours of one cell, you need to disturb the others. 

# plot by doublet score
plotUMAP(sce, colour_by="scDblFinder.score") +
  plotTSNE(sce, colour_by="scDblFinder.score") +
  plot_layout(guides = "collect")

# filter out bad cells
sce <- sce[,sce$scDblFinder.class!="doublet" & !sce$qc.out]

Both representations are capturing the same time of relationships (you can’t not really observe it here but yes in the plot he presented with more colours)

Batch correction

First we want to check whether there’s a need for integration/batch correction. This can be assessed by the extent to which the samples are well-mixed.

# check mixing:
plotTSNE(sce, colour_by="group_id") +
  plotTSNE(sce, colour_by="sample_id") +
  plot_layout(guides = "collect")

## a much better readout (but long to compute) for cell mixing would be given 
## by the CellMixS package:
# sce <- cms(sce, k=50, group = "sample_id", BPPARAM=MulticoreParam(6))
# plotTSNE(sce, colour_by="cms_smooth")
# batch corrected using the mutual nearest neighbors
sce2 <- fastMNN(sce, batch=sce$sample_id, BNPARAM=AnnoyParam())
# (see Harmony or Seurat integration for "stronger" integration methods)
# we take the corrected PCA
reducedDim(sce, type="PCA") <- reducedDim(sce2)[,1:20]
# and recompute the 2d projections
sce <- runTSNE(sce, dimred="PCA", BNPARAM=AnnoyParam())
sce <- runUMAP(sce, dimred="PCA", n_neighbors=25, BNPARAM=AnnoyParam())

plotUMAP(sce, colour_by="group_id") + 
  plotUMAP(sce, colour_by="sample_id")

Integration from the SEURAT packaged or harmony are useful to stronger batch-effect corrections.

Clustering

# BiocNeighbors
g <- buildKNNGraph(sce, BNPARAM=AnnoyParam(), use.dimred="PCA", k=30)
sce$cluster <- as.factor(cluster_leiden(g, objective_function="modularity",
                                        n_iterations=20)$membership)
table(sce$cluster)

   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
 176  565  272  625  898  399 1252 1597  301 1043  955 1573  380   94   49  591 
  17 
 689 
# we could play with the `resolution_parameter` of `cluster_leiden` 
# to decide on the granularity

plotTSNE(sce, colour_by="cluster", text_by="cluster") +
  plotUMAP(sce, colour_by="cluster", text_by="cluster")

On the slide 10/11, we can see that as increasing resolution, some clusters divide into others. That indicates that clustering might not have the best quaility.

Cluster abundances across samples

ca <- table(cluster=sce$cluster, sample=sce$sample_id)
ggplot(as.data.frame(ca), aes(reorder(sample), cluster, fill=Freq)) +
  geom_tile() + geom_text(aes(label=Freq)) +
  scale_fill_viridis_c()

In the case of our data, we don’t really expect a lot of differencies in abundance.

Cluster annotation

De-novo marker identification

# we identify genes that differ between clusters:
mm <- scran::findMarkers(sce, groups=sce$cluster, test.type="binom",
                         BPPARAM=MulticoreParam(6))
# we select the top 5 markers by cluster:
markers <- unique(unlist(lapply(mm, FUN=function(x){
  head(row.names(x[x$FDR<0.01,]),5)
})))

Known markers

This comes from literature information. What is a good marker as protein and stains might not be a good marker at RNA level.

genes <- list(
  astrocytes = c("Aqp4", "Gfap", "Fgfr3","Dio2"),
  endothelial = c("Cldn5","Nostrin","Flt1"),
  microglia = c("C1qb","Tyrobp","P2ry12", "Csf1r"),
  neuronal = c("Snap25", "Stmn2", "Syn1", "Rbfox3"),
  excNeuron = c("Slc17a7","Camk2a","Grin2b","Fezf2"),
  inhNeuron = c("Gad1","Lhx6","Adarb2"),
  oligodendro = c("Opalin","Plp1","Mag","Mog"),
  OPC = c("Pdgfra","Sox6","Bcan")
)

# since the row.names of the object have also the ensembl id, we find the matching row names for each gene:
km <- lapply(genes, FUN=function(g) grep(paste0(g, "$", collapse="|"), rownames(sce), value=TRUE))

Pseudo-bulk aggregation

# muscat
# mean logcounts by cluster:
pb <- aggregateData(sce, "logcounts", by=c("cluster"), fun="mean")
assayNames(pb) <- "logcounts"
assays(pb)$propOfMax <- exp(logcounts(pb))/rowMaxs(exp(logcounts(pb)))
rowData(pb)$marker4 <- NA
rowData(pb)[unlist(km),"marker4"] <- rep(names(km),lengths(km))

sechm(pb, c(unlist(km)), assayName="logcounts", gaps_row="marker4",
      show_colnames=TRUE, do.scale=TRUE, breaks=1, row_title_rot=0) + 
  sechm(pb, c(unlist(km)), assayName = "propOfMax", show_colnames=TRUE,
        do.scale=FALSE, hmcols=viridis::viridis(100), 
        row_names_gp=gpar(fontsize=9))

# heatmap for the de-novo markers:
sechm(pb, markers, assayName = "propOfMax", show_colnames=TRUE,
        do.scale=FALSE, hmcols=viridis::viridis(100), 
        row_names_gp=gpar(fontsize=9))

For simplicity, here rather than working on all clusters we’ll group them by broad cell class. As a very rough approximation to manual annotation, we will assign clusters to the cell type whose markers are the most expressed:

# we get rid of the unspecific neuronal markers:
km2 <- km[names(km)!="neuronal"]
# we extract the pseudo-bulk counts of the markers for each cluster
mat <- assay(pb)[unlist(km2),]
# we aggregate across markers of the same type
mat <- aggregate(t(scale(t(mat))),
                 by=list(type=rep(names(km2), lengths(km2))),
                 FUN=sum)
# for each column (cluster), we select the row (cell type) which has the maximum aggregated value
cl2 <- mat[,1][apply(mat[,-1], 2, FUN=which.max)]

# we convert the cells' cluster labels to cell type labels:
sce$cluster2 <- cl2[sce$cluster]
table(sce$cluster, sce$cluster2)
    
     astrocytes endothelial excNeuron inhNeuron microglia oligodendro  OPC
  1           0           0       176         0         0           0    0
  2           0           0       565         0         0           0    0
  3           0           0         0         0       272           0    0
  4           0           0       625         0         0           0    0
  5           0           0       898         0         0           0    0
  6           0         399         0         0         0           0    0
  7           0           0      1252         0         0           0    0
  8        1597           0         0         0         0           0    0
  9           0           0         0       301         0           0    0
  10          0           0         0      1043         0           0    0
  11          0           0       955         0         0           0    0
  12          0           0         0         0         0        1573    0
  13          0         380         0         0         0           0    0
  14          0          94         0         0         0           0    0
  15          0           0         0         0         0          49    0
  16          0           0         0       591         0           0    0
  17          0           0         0         0         0           0  689
plotUMAP(sce, colour_by="cluster2", text_by="cluster2")

# we aggregate again to pseudo-bulk using the new clusters
pb <- aggregateData(sce, "logcounts", by=c("cluster2"), fun="mean")
assayNames(pb) <- "logcounts"
# and we plot again the expression of the markers as a sanity check
rowData(pb)$marker4 <- NA
rowData(pb)[unlist(km),"marker4"] <- rep(names(km),lengths(km))
sechm(pb, c(unlist(km)), assayName="logcounts", gaps_row="marker4",
      show_colnames=TRUE, do.scale=TRUE, breaks=1, row_title_rot=0)

plotTSNE(sce, colour_by="cluster2", text_by="cluster2")

saveRDS(sce, file="week13.SCE.processed.rds")

Differential state analysis

# we aggregate by cluster x sample to perform pseudo-bulk differential state analysis
sce <- muscat::prepSCE(sce, kid="cluster2", sid = "sample_id",
                       gid = "group_id")
pb <- aggregateData(sce)

pbMDS(pb)

# we run edgeR on each cluster and extract the results
res <- pbDS(pb)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |======================================================================| 100%
res2 <- dplyr::bind_rows(res$table[[1]])

# volcano plots
ggplot(res2, aes(logFC, -log10(p_adj.loc), colour=p_adj.loc<0.05)) + geom_point() + facet_wrap(~cluster_id)

# top genes in a given cell type
pbHeatmap(sce, res, k="astrocytes", top_n = 5)

# we extract all differentially-expressed genes:
degs <- unique(res2[res2$p_adj.loc<0.05,"gene"])

# we flatten the pb object (putting all cell types in the same assay) and calculate foldchanges
pb2 <- pbFlatten(pb)
# we add a logFC assay:
pb2 <- sechm::log2FC(pb2, fromAssay="logcpm", controls=pb2$group_id=="WT", by=pb2$cluster_id)
# we reorder
pb2 <- pb2[,order(pb2$cluster_id, pb2$group_id)]
# we plot a heatmap of the logFC of the top 200 genes across all cell types
sechm(pb2, head(degs,200), assayName="log2FC", gaps_at="cluster_id",
      column_title_gp=gpar(fontsize=9))

The number of differentially expressed genes it is a bad signature for response to treatment when we have substantial differencies in abundances.

Session info

sessionInfo()
R version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Zurich
tzcode source: internal

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] patchwork_1.1.3             sechm_1.9.5                
 [3] ComplexHeatmap_2.16.0       igraph_1.6.0               
 [5] BiocNeighbors_1.18.0        BiocParallel_1.34.2        
 [7] muscat_1.14.0               scDblFinder_1.14.0         
 [9] batchelor_1.16.0            scater_1.28.0              
[11] ggplot2_3.4.4               scran_1.28.2               
[13] scuttle_1.10.3              SingleCellExperiment_1.22.0
[15] SummarizedExperiment_1.30.2 Biobase_2.60.0             
[17] GenomicRanges_1.52.1        GenomeInfoDb_1.36.4        
[19] IRanges_2.34.1              S4Vectors_0.38.2           
[21] BiocGenerics_0.46.0         MatrixGenerics_1.12.3      
[23] matrixStats_1.2.0          

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.21          splines_4.3.1            
  [3] BiocIO_1.10.0             bitops_1.0-7             
  [5] tibble_3.2.1              XML_3.99-0.16            
  [7] lifecycle_1.0.4           Rdpack_2.6               
  [9] edgeR_3.42.4              doParallel_1.0.17        
 [11] globals_0.16.2            lattice_0.22-5           
 [13] MASS_7.3-60               backports_1.4.1          
 [15] magrittr_2.0.3            limma_3.56.2             
 [17] rmarkdown_2.25            yaml_2.3.8               
 [19] metapod_1.8.0             sctransform_0.4.1        
 [21] cowplot_1.1.2             minqa_1.2.6              
 [23] RColorBrewer_1.1-3        ResidualMatrix_1.10.0    
 [25] abind_1.4-5               zlibbioc_1.46.0          
 [27] Rtsne_0.17                EnvStats_2.8.1           
 [29] glmmTMB_1.1.8             purrr_1.0.2              
 [31] RCurl_1.98-1.13           seriation_1.5.4          
 [33] circlize_0.4.15           GenomeInfoDbData_1.2.10  
 [35] ggrepel_0.9.4             pbkrtest_0.5.2           
 [37] irlba_2.3.5.1             listenv_0.9.0            
 [39] dqrng_0.3.2               parallelly_1.36.0        
 [41] DelayedMatrixStats_1.22.6 codetools_0.2-19         
 [43] DelayedArray_0.26.7       tidyselect_1.2.0         
 [45] shape_1.4.6               farver_2.1.1             
 [47] lme4_1.1-35.1             ScaledMatrix_1.8.1       
 [49] viridis_0.6.4             TSP_1.2-4                
 [51] GenomicAlignments_1.36.0  jsonlite_1.8.8           
 [53] GetoptLong_1.0.5          randomcoloR_1.1.0.1      
 [55] iterators_1.0.14          foreach_1.5.2            
 [57] tools_4.3.1               progress_1.2.3           
 [59] Rcpp_1.0.11               blme_1.0-5               
 [61] glue_1.6.2                gridExtra_2.3            
 [63] xfun_0.41                 mgcv_1.9-0               
 [65] DESeq2_1.40.2             ca_0.71.1                
 [67] dplyr_1.1.4               withr_2.5.2              
 [69] numDeriv_2016.8-1.1       fastmap_1.1.1            
 [71] boot_1.3-28.1             bluster_1.10.0           
 [73] fansi_1.0.6               caTools_1.18.2           
 [75] digest_0.6.33             rsvd_1.0.5               
 [77] R6_2.5.1                  colorspace_2.1-0         
 [79] Cairo_1.6-2               gtools_3.9.5             
 [81] RhpcBLASctl_0.23-42       utf8_1.2.4               
 [83] tidyr_1.3.0               generics_0.1.3           
 [85] variancePartition_1.30.2  data.table_1.14.10       
 [87] rtracklayer_1.60.1        prettyunits_1.2.0        
 [89] htmlwidgets_1.6.4         S4Arrays_1.0.6           
 [91] uwot_0.1.16               pkgconfig_2.0.3          
 [93] gtable_0.3.4              registry_0.5-1           
 [95] XVector_0.40.0            remaCor_0.0.16           
 [97] htmltools_0.5.7           TMB_1.9.10               
 [99] clue_0.3-65               scales_1.3.0             
[101] png_0.1-8                 knitr_1.45               
[103] rstudioapi_0.15.0         reshape2_1.4.4           
[105] rjson_0.2.21              curl_5.2.0               
[107] nlme_3.1-164              nloptr_2.0.3             
[109] GlobalOptions_0.1.2       stringr_1.5.1            
[111] KernSmooth_2.23-22        parallel_4.3.1           
[113] vipor_0.4.5               restfulr_0.0.15          
[115] pillar_1.9.0              vctrs_0.6.5              
[117] gplots_3.1.3              BiocSingular_1.16.0      
[119] beachmat_2.16.0           cluster_2.1.6            
[121] beeswarm_0.4.0            evaluate_0.23            
[123] isoband_0.2.7             mvtnorm_1.2-4            
[125] cli_3.6.2                 locfit_1.5-9.8           
[127] compiler_4.3.1            Rsamtools_2.16.0         
[129] rlang_1.1.2               crayon_1.5.2             
[131] future.apply_1.11.0       labeling_0.4.3           
[133] plyr_1.8.9                ggbeeswarm_0.7.2         
[135] stringi_1.8.3             viridisLite_0.4.2        
[137] lmerTest_3.1-3            munsell_0.5.0            
[139] Biostrings_2.68.1         aod_1.3.3                
[141] V8_4.4.1                  Matrix_1.6-4             
[143] hms_1.1.3                 sparseMatrixStats_1.12.2 
[145] future_1.33.0             statmod_1.5.0            
[147] rbibutils_2.2.16          broom_1.0.5              
[149] xgboost_1.7.6.1